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Abstract 

Experimental measurements of properties of the large-scale circulation (LSC) in turbulent con- 
vection of a fluid heated from below in a cylindrical container of aspect ratio one are presented 
and used to test a model of diffusion in a potential well for the LSC. The model consists of a pair 
of stochastic ordinary differential equations motivated by the Navier-Stokes equations. The two 
coupled equations are for the azimuthal orientation 9q^ and for the azimuthal temperature ampli- 
tude 5 at the horizontal midplane. The dynamics is due to the driving by Gaussian distributed 
white noise that is introduced to represent the action of the small-scale turbulent fluctuations on 
the large-scale flow. Measurements of the diffusivities that determine the noise intensities are re- 
ported. Two time scales predicted by the model are found to be within a factor of two or so of 
corresponding experimental measurements. A scaling relationship predicted by the model between 
6 and the Reynolds number is confirmed by measurements over a large experimental parameter 
range. The Gaussian peaks of probability distributions p{5) and p{9q) are accurately described by 
the model; however the non-Gaussian tails of p{5) are not. The frequency, angular change, and 
amplitude bahavior during cessations are accurately described by the model when the tails of the 
probability distribution of 5 are used as experimental input. 
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Chicago, IL 60637 



I. INTRODUCTION 



Turbulent convection is one of the outstanding unsolved problems of classical physics 
(for reviews, see for example Refs. l|, [2I, S]). The problem of Rayleigh-Benard convection 
(RBC) consists of a fluid sample heated from below. The heat input causes the emission of 
volumes of hot fluid known as "plumes" from a bottom thermal boundary layer that rise due 
to buoyancy, while cold plumes emitted from a top boundary layer sink. The experiments 
are done in cylindrical containers with an aspect ratio T = D / L \ (Lis the height and D 
is the diameter of the sample). In the turbulent regime of F = 1 samples, the plumes drive 
a large-scale circulation (LSC), also known as the "mean wind", which is oriented nearly 



vertically with up-flow and down-flow on opposite sides of the sample 
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The wind in turn carries the plumes, primarily up one side and down the other. The 



dy„a,„,cs of the LSC ...c.ude oscillations HttauHflnnHQ^in wh.ch^e 

orientation of the upper half of the LSC oscillates out of phase with the lower half [9|, |20J. 
The LSC breaks the rotational symmetry of a cylindrical sample and its circulation plane 
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must somehow choose an azimutha 
spontaneous diffusive meandering 
both by azimuthal rotations jl^ . I23I ]. and 
and restarts in a random new orientation 



orientation. This orientation has been found to undergo 
2^ 
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2J]. The LSC also undergoes re-orientations 



by: 
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cessations in which the LSC slows to a stop 



251]. On longer time scales. Earth's Coriolis 



force was found to cause a net rotation of the LSC orientation on average once every 3 days, 
and to align the LSC in a preferred orientation close to West [2^. These LSC dynamics 
observed in experiments may be related to some natural convection dynamics. For example, 
it is possible that cessations of the flow in the Earth's outer core are responsible for changes 
in the orientation of Earth's magnetic fleld 26|]. Reversals are known to occur in the wind 



direction in the atmosphere 27|]. Torsional oscillations are observed in the solar convection 
zone 28| . 

Two stochastic models of LSC flow-reversal have been proposed in the literature 29|, l30 |. 
They treated diffusion of the LSC strength in a potential well, but there was no physical 
motivation for the shape of the potential that was used (which differs from ours) and the 
model parameters were chosen phenomenologically. No azimuthal degree of freedom was 



included, and thus only genuine reversals of the LSC 



which are now known to be very rare 



events) could be produced. Two other models [3l|, l32] describe the LSC with deterministic 
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differential equations that have chaotic solutions. They have their roots in the Navier-Stokes 



equations and retain terms that are argued to be physically important. One of them [31 1 
again is lacking the azimuthal degree of freedom that is so important to the LSC dynamics 
found in experiment. The other [32] is based on an exact solution of the Boussinesq equations 
in the inviscid and unforced limit, but employs physically unrealistic boundary conditions 
and adds dissipation a posteriori. It lacks the stochastic character found in experiment, and 
requires the arbitrary adjustment of parameters to yield the chaotic solutions that might be 
compared with observations. 

In order to improve upon the state of the field described above, we presented briefly in 
Ref. [33!] a stochastic model of the LSC that was motivated by the physically relevant terms 
of the Navier-Stokes (NS) equations for RBC It is in the same spirit as a model for the effects 



of Earth's Coriolis force on the flow 



24l |. and we will show in a subsequent paper that the 



Coriolis-force model is consistent with the strong-damping limit of the current model. The 
model consists of two coupled stochastic ordinary differential equations (ODEs): one for the 
strength of the LSC represented by an amplitude 6 of the azimuthal temperature variation 
at the horizontal mid-plane of the sample, and the other for the azimuthal LSC orientation 
^0- For 6 it leads to diffusive motion in a potential which has a minimum at 6q > and a 
maximum at 5 = 0. On the rare occasions when the diffusion reaches (or comes close to) the 
maximum, then a cessation has occurred. The shape of that potential follows from taking a 
volume average of contributions from buoyancy and from the viscous drag on the walls. The 
equation for 6q contains a nonlinear coupling to the 5-equation which is proportional to S, 
which comes from the advective term in the NS equation, and which represents the angular 
momentum of the circulation. Thus, when 6 is large, the angular momentum is large and the 
LSC orientation is relatively immune to re-orientation by the stochastic forces. On the other 
hand, near cessations where 6 becomes small, re-orientations are relatively easily achieved 
by the fluctuating background. The stochastic driving represents phenomenologically the 
action of the small-scale turbulent fluctuations on the large-scale circulation. The strength 
of these fluctuating forces is obtained from experiment by measuring the diffusivities of the 
relevant variables. The model was solved numerically, and the numerical solutions were 
found to reproduce quite well the cessations and rotations, as well as the diffusive azimuthal 



meandering, that had been observed in experiments 



The model also yielded 



probability distributions for the angular displacement during cessations and rotations that 
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were in quite good agreement with the measurements |23l. |25|. 

In the present paper we present a more detailed derivation of and motivation for the 
model. Then we derive new analytic results by a lowest-order expansion of the potential for 
the 6 equation about its minimum, and by linearizing the coupling term of the equation. 
We present experimental measurements of the model parameters, as well as of the Rayleigh- 
number dependences of several quantities, and for the most part find reasonable agreement of 
these results with the model predictions. A notable exception is the shape of the 5-potential 
near 6 = 0. We suggest that a possible explanation for this disagreement may be found in 
the neglect by the model of diffusive heat transport across the top and bottom boundary 
layers which can be expected to become significant when the LSC amplitude becomes small. 

In the next section we discuss the experiment. This is kept brief because much of it 
was done before 
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34j ]. In Sect. IIII Al we give a detailed derivation of and motivation 
for the model. A linearized version of the model is derived in Sect. IIII B[ The potential 
for the 6 equation, and its linearized version, are discussed in Sect. IIII CI In Sect. HVl we 
present new experimental measurements of the model parameters. First measurements of 
the Reynolds numbers are discussed. Then the diffusivities corresponding to stochastic 
fluctuations are presented. Next, probability distributions, power spectra, and correlations 
of 6 and 6q are given. This is followed in Sect |V] by a comparison between model predictions 
for various Rayleigh-number dependences and measurements of the model time scales, of 
6o, of the diffusivities, and of further interesting dimensionless parameters. In Sect. I VI I we 
use the experimental measurement of the tails of the probability distribution of 6 as an 
experimental input to the model to predict the behavior of cessations more accurately. We 
predict the duration, frequency, and net angular change during cessations and compare these 



values to the experimental results of Ref. 



23|. In Sect. IVIII we suggest a modification of 



the model that would add a third equation and that would be expected to account for the 
small-5 behavior without the empirical experimental input. 

II. THE EXPERIMENT 

In order to measure the model parameters, experimental data from two cylindrical sam- 
ples with aspect ratio F 1 were used. These were the medium and large sample described 



in detail elsewhere 
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34| . The samples had copper top and bottom plates and a plexiglas 
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side wall. The medium sample had a diameter D = 24.81 cm and a height L = 24.76 cm, 
and the large sample had D = 49.67 cm and L = 50.61 cm. The Rayleigh number R is given 
by 

R^^i^ (1) 

where a is the isobaric thermal expansion coefficient, g the acceleration of gravity, AT the 
applied temperature difference, k the thermal diffusivity, and u the kinematic viscosity. By 
applying a temperature difference 0.5 K < AT < 20 K between the bottom and top plates, 
with two samples of different heights L, a Rayleigh number range of 3 x 10^ < R < 10^^ 
could be covered. The Prandtl number a is given by 

(2) 

Each sample was filled with water and the average temperature between the bottom and top 
plates usually was kept at 40.0° C, giving cr = 4.38. Some measurements were made at other 
temperatures and permitted a change of a over the range 3.3 ^ cr ^ 5.5. Measurements were 
made with thermistors placed into blind holes drilled into the side wall from the outside so 
they did not interfere with the flow. There were eight thermistors at the mid-height of the 
side wall, equally spaced azimuthally. The LSC carried warm fluid from the bottom plate 
up one side of the sample, which cooled when it passed the top plate and went down on the 
opposite side of the sample. The temperature proffie 

T = To + (5cos(0o-^) (3) 
was fit to the temperature measurements where 5 characterizes the strength of the LSC and 



^0 is the orientation of the LSC 23|, |35||. Fits of the temperature measurements every 2.5 



seconds provide time series of b and^. These time series contained the diffusive dynamics 



of the LSC and the re-orientations 



23 



25| . The model parameters can be extracted from 



the time series as described in Sect. [IVl 
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III. THE MODEL 



Derivation 



The model presented in Ref. 



33| is reproduced here in more detail. The LSC strength 



can be described by the velocity component u,^. Here (j) is an angle in the vertical circulation 
plane of the LSC as shown in Fig. [H and describes the flow in the absence of azimuthal 
motion. One expects the acceleration to be due to a balance between buoyancy and drag 
forces. The pressure term primarily provides the inward force to keep the LSC in a loop, 
and it is not expected to contribute a net force in the 0-direction. Thus only the buoyancy 
and viscous drag terms are included on the right hand side of the NS equation for m^, and 
we neglect the nonlinear term: 



ii^ = ga(T - Tq) + uV^u^ . 



(4) 



To obtain a model in the form of an ODE that describes the flow with only a few variables, 



a global average is taken over the field variables that retains the essential physics o: 
LSC. This average can be carried out using the experimental observation 2 
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the 



36| 



that the temperature of the LSC at the side wall at mid-height is given by Eq. [31 The profile 
is interpolated to be 



2rS 

T{r,e) = To + —cos{do 



(5) 



where r is the radius measured from the cylinder axis. An approximately linear radial 
variation was found experimentally in Ref. Sj. The buoyancy acts on the entire LSC. It 
enhances the LSC in proportion to its vertical component. Thus, to approximately ac- 
count for this, we multiply the volume average by a factor of 1/2. The buoyancy term can 
now be approximated by an integral over the container volume V using Eq. [5] to obtain 
l/(2\^) Jy ga(T — Tq)S{9 — 9Q)dV = 2ga6/{3'K) where 5* is a step function with S* = 1 for 
1^0 — ^1 < 7r/2 and S = —1 otherwise. Note that when the volume averages are taken, 
assumptions about the geometry of the flow only affect coefficients by factors of order one, 
and do not influence the functional form of the results. 
It was argued in Ref. 



23| based on the experiments of Ref. [lOj that the azimuthal velocity 



6 



profile is close to a step function, so the bulk velocity profile is assumed to be 

u,ir,e) = ^ (6) 

where U is the maximum vertical speed near the side wall. Here again the linear variation 
with r is supported by experimental results [s, [igI. The drag is assumed to occur in the 
viscous boundary layers, where V^u^ ~ —U/X^ (A is the viscous boundary-layer width). 
The viscous boundary layers on the side wall and plates occupy a fraction of the container 
equal to 6X/L, so the volume average of the drag is —QuU/{XL). The viscous boundary- 

— 1/2 

layer width is assumed to follow the Prandtl-Blasius form A = LR^/ /2 with a fluctuating 
Reynolds number R^^i = UL/u. Although this must be regarded as an approximation, the 
Prandtl-Blasius form for the boundary layer has worked remarkably well in previous models 



(for example, 24] )• 1^ ^Iso has been very successful in predicting the dependence of the 



Reynolds number on the Rayleigh number [37[, and in treating non-Boussinesq effects on 



391] . With this form the damping term 



the Nusselt number and the center temperature [38|, 
becomes nonlinear in U. A volume average of the acceleration term using Eq. [6] results in 
(l/V) Jyii^dV = 2U /3. Combining these results gives the volume-averaged equation 



2U_ _ 2ga5_ _ 12i^^/^U^/^ 

This equation has two variables, 6 and U, but we only have experimental measurements of 6. 
In order to compare the model to current data, it is assumed that the temperature amplitude 
6 is instantaneously proportional to the speed U, since both variables are measures of the 
LSC strength. This assumption is consistent with simultaneous velocity and temperature 
measurements at the same point at the mid-height near the side wall which gave a correlation 
of 0.8 [4^. Two-dimensional direct numerical simulations also found that on average 6 is 
proportional to U, although instantaneously both fluctuate, with the same correlation of 0.8 
4l| . The proportionality constant relating 6 to U must satisfy the time-average 

2ga{6) Uv^Rll^ 

of Eq. [71 This fixes the proportionality at 

2ga5 _ UvURy^ 
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Note that Eq. [8] forces the sum of the powers of U and Re to be 3/2, and so the assumption 
fixes 5 (X URy^. Equation [9] is substituted into Eq. [7] and a stochastic term fs{t) is added 
to represent the influence of the small-scale turbulent fluctuations on the large-scale flow to 
get the Langevin equation 

5= jr + ^(^) • (10) 

Using Eqs. [1] and [2] one finds the constant 

-^0 = ^ (11) 

and the time scale 



The stochastic term is assumed to be Gaussian distributed white noise with zero mean. 
These properties will be explored in detail in Sect. HVl 

Equation has two fixed points, one unstable at 5 = and one stable when 6 = 6o- 
Thus, in the absence of fluctuations, 60 can be interpreted as the steady-state amplitude. 
The stochastic equation reproduces some of the important behavior of the LSC When a 
temperature difference is applied to generate buoyancy, the LSC will start to grow due to 
the instability at 5 = until it reaches the stable fixed point at 6 = 6q. If the fluctuations 
are small, 6 spends most of its time meandering near the stable fixed point at 6q, and if the 
fluctuations are large enough the LSC occasionally undergoes a cessation when fluctuations 
drive 6 close to 0. 

Other models of the LSC dynamics 0, [s^ have used an equation for the LSC strength 
similar to Eq. [TOl but assumed an exponent of 3 for the damping term instead of 3/2. The 
essential physics of Eq. [10] is in the (in)stability behavior of the fixed points. Equations like 
Eq. [To] have one unstable and one stable fixed point as long as the damping exponent is 
greater than one. Thus in a qualitative sense Eq. [10] has a behavior similar to that of the 
earlier models. Since 6 is chosen to be non-negative and reversals are accounted for by a 
change in orientation^, there is no need to restrict the exponents to odd integers as in the 



^ The parameters of the temperature profile T = Tq + S cos{9 — 60) are not uniquely determined, because the 
change 6 ~6 is equivalent to Oq ^ 9o± n. Thus S is chosen to be always non-negative for uniqueness. 
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other models. Since the 3/2 power came from our choice of drag law and scaling relationship 
between 5 and f/, the phenomenology of the model is also robust to these choices as long 
and the exponent of the drag term is greater than that of the buoyancy term. 

The second Langevin equation describes the azimuthal motion of the LSC The only 
driving force for azimuthal motion in a symmetric system is turbulent fluctuations, and 
damping can come from either viscosity or rotational inertia. Thus only the drag term is 
kept on the right-hand-side of the NS equation in the azimuthal coordinate: 

iie + u - Vu0 = uV'^ue ■ (13) 

Again, the components of the nonlinear term are negligible except the one corresponding to 
the rotational inertia of the LSC in the ^-coordinate [u- V)ue ~ {u^/r)dug/ dcj) ~ U9q. This 
can be physically understood in terms of the dynamics of a rigid rotator: a rotating body has 
some stability due to its angular momentum that resists a torque in an orthogonal direction. 
The torque in the ^-direction is equal to I9q = Lg where Lg and are the angular momenta 
in the respective coordinates, and I is the moment of inertia. Since the LSC nearly fills the 
container, I is assumed to be the same around both axes of rotation. For a differential torque 
applied to a rigid rotator, a change in orientation is more difficult when there is rotation 
in a perpendicular direction d9o = dLg/L^, or Lg = L^Oq, where = lu^/r ^ 2UI/L. 
Combining these equations yields the inertial contribution to acceleration 6q = 2U6q/L. 
Using the approximation that viscous drag occurs mainly in the boundary layers, one has 
{h'V'^ug)v ~ v{L6q/2)/\^ X 6A/L. The volume average of Eq. [T3l is 

3 3 L ^ ^ 

— 1/2 

The ratio of the viscous drag term to the angular-momentum damping is equal to 9R^ ^ . 
At i?e,i = 3700 (the Reynolds number at i? = 1.1 x 10^°) for example, this yields ~ 0.15 for 
this ratio. Since rotational inertia damps the azimuthal motion much more than the viscous 
drag across the boundary layer near the side wall, the viscous damping term of the azimuthal 
equation is neglected from now on. The azimuthal speed is generally small compared to the 
LSC speed so the effect of rotational inertia is much larger on the azimuthal coordinate than 
on the LSC strength, which is why the nonlinear term could be ignored in Eq. HI Converting 
the remaining terms from U to 6 using Eq. [9|, and adding another stochastic term fg{t) 
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representing turbulent fluctuations, gives 



with the constant time scale 



Oo = -^ + feit) (15) 



(16) 



" 2uRe 

The two stochastic ODEs Eqs. [10] and [15] compose the model for the LSC dynamics. The 
stochastic terms fs(t) and fg(t) that drive the dynamics of the system are presumed to 
originate from the small-scale turbulent background fluctuations. We made some extremely 
simplifying assumptions about them because they cannot be isolated from the dynamical 
system to independently measure their properties. They are assumed to be Gaussian dis- 
tributed, uncorrelated white noise. Then the diffusivities Ds and Dg are the only parameters 
required to describe them, and these can be estimated from experimental data. The ade- 
quacy of these assumptions is tested in the context of the model in Sect. IIVI There are 
direct model predictions for 6q and the time scales ts and Tq. Methods for experimentally 
obtaining all of the parameter values and testing the model will be covered in Sect. llVi 



B. Linear approximation 

When fluctuations of 6 are small, or 5 ^ 6n, a linear approximation to the Langevin 



equations can be made. For the data of Ref. [23|] this approximation is satisfied most of 
the time, but not during the quite rare cessations. Thus the linearized version of the model 
should apply to long-term averages of data, but the full non-linear model must be used to 
study cessations. 

Starting with Eq. [TO] for 6, we expanded around the stable fixed-point solution 5 = 5o by 
rewriting 6 = 6o + e. Assuming e <^ 6o and thus keeping up to the first order term in the 
expansion one has (5o + e)3/2 ^ ^3/2 ^ 3g^/2. Using 6{6o) = 0, Eq. [10] simplifies to a linear 
equation 

e = -^ + fs{t). (17) 
Its 

Equation [15] can also be linearized near the stable fixed point by setting 6 = 6o to obtain 
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Oo = -- + fe{t). (18) 

The typical size of fluctuations in S can be measured by the ratio of the variance of 6 to the 
average of 6 squared, or crf/^Q ^ 0.07 for the large sample data and smaller in the medium 
sample (see Sect. IIVCI) . Since this ratio is small, the effect of the variable nature of the 
damping is on average small. 



C. The potential wells 

It is useful to think of Eq. [TO]in terms of diffusion in a potential well, as in the Arrhenius- 



Kramers problem 



42l | . The potential is defined by = — / 6dd6 where 5^ is the deterministic 



part of Eq. [TOl Thus 

This potential is shown in Fig. [21 Its minimum is at 5 = 5o, and due to the stochastic 
term fs{t) the value of 6 fluctuates around 6q in the bottom of the well. A cessation occurs 
when the LSC amplitude drops to zero, or when fiuctutations in 6 cross the potential barrier 
AV^ = Vs{0) — Vs{So). The linearization of the potential well near the stable fixed point 
gives 

K = Vsi6o) - I idde = Vsi6o) + , (20) 

where q is the deterministic part of Eq. [T71 This potential is also shown in Fig. [2l The 
minima of both potentials overlap and have the same curvature, so the dynamics for small 
fluctuations will be the same for both; but the non-linear potential Vs is skewed to bias 
fluctuations towards small 6 relative to the parabolic potential K- 

Similarly, the dynamics of given by Eq. [TSl can be thought of as diffusion in the 
potential well 

Ve = -1 9od9o = elS/{2SoT^) . (21) 

In this case the well is parabolic, with a minimum at ^ = 0. Thus the mean of 6 will be zero, 
with fluctuations symmetrically about this value. However, the well curvature varies with 6. 
Thus, reorientations occur when fluctuations take 60 far away from the average of zero, and 
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this is more likely when S, and thus the curvature of the potential, are small. The extreme 
case is a cessation, which occurs when 6 is close to zero and the potential Vg is nearly flat. 



IV. MEASURING THE MODEL PARAMETERS 



A. The Reynolds number and the mean temperature ampHtude 6q 

The model of the large-scale circulation requires several input parameters. The Reynolds 

43] ) and 



number i?e was measured in many experiments (for a recent summary, see Ref. 
desribed by theoretical models [sT^ . Values are taken from temperature-correlation functions 
corresponding to the plume turnover that came from the same apparatus as the current mea- 



surements and were reported in Ref. [43j, and will not be distinguished from different mea- 



sures of the Reynolds number that were discussed in Ref. 43| because these differences are 
small compared to the needs of an order-of-magnitude model. The flxed-point temperature- 
amplitude ^0 can be approximated by an average over a time series: 6o ~ (6). Due to the 
asymmetry of Eq. [19] around 6 = 6q, (6) is slightly less than 6q, but the difference is a small 
fraction of 6o and so will be ignored. 



B. The diffusivities and time scales 

For diffusive fluctuations, the mean-square change {{dx)"^) of a variable x over a time 
interval dt is a linear function of dt, and the diffusivity is deflned by the equation 
{{dxY) = Dj-dt. For stepwise numerical simulations, the stochastic terms fx{t) in the model 
have a variance D^/h where h is the time step of the numerical integration. 



1. Diffusion of the amplitude 6 

A plot of {{d5Y) = {[6(t + dt) — S{t)]'^) as a function of the time interval dt is shown in 
Fig. [3] for i? = 1.1 X 10^° as an example (data at numerous other values of R behave in a 
similar manner). The equation {{dS)"^) = Dsdt was flt to the data for 30s < dt < 80 s to 
obtain Ds = 6.4 ■ 10~^K'^/s. Although this linear relationship is characteristic of diffusion, 
the variance of the change in 6 saturates at a constant value for large dt . This happens 
because the diffusion of 6 occurs in the potential well given by Eq. \TU\ and so is bounded 
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to a finite range. For the linearized model Eq. [T71 the potential well is parabolic so p{S) is 
Gaussian with variance cr| = tsDs (see Sect II V CI for a derivation). The long-term variance of 
the change in 6 is given by \imdt^^{[6{dt) -6{0)]^) = \imdt^^{[S{dt)Y -2S{dt)6{0) + [6{0)]'^). 
Since 6{dt) and 5(0) are uncorrelated for large dt, but each have variance erf, it follows for 
large dt that {[d5{dt)]^) = 2oj = 2tsDs. Since Dg is measured from the slope of {[dS^dt)]"^) 
for intermediate dt in Fig. [3], one can obtain ts from the ratio of measured values to get 
Ts = 47 s. The transition between the two scaling regimes occurs at dt = 2ts. Thus 2ts is the 
time scale over which the amplitude retains some correlation. The parameters determined 
here and below for i? = 1.1 x 10^° are summarized in Table [B The dependences of the 
measured diffusivities and time scales on R are shown in Sect. El 

2. Diffusion of 9q 

Figure m shows the mean-square change in azimuthal rotation rate {{dOY) = {6o{t + dt) — 
0o{t)) as a function of the time interval dt. Here 6o{t) = [6o{t + dt/2) — 6o{t — dt /2)]/ dt. Also 
plotted is the same quantity derived from measurements of ^o(^) that are restricted to the 
range 0.95o < S < I.ISq near the stable fixed point. The equation {d6l) = Dgdt is fit to the 
latter data to obtain Dg = 2.9 x 10^^ rad^/s'^ for i? = 1.1 x 10^°. The diffusivity is calculated 
from the data with 6 close to 6o so that it can be analyzed according to the linear prediction 
Eq. [iHl The difference between the two results is small in any case. The plot of {[dOoldt)]"^) 
has a plateau because the damping term causes 6q to be bounded. Since the linear azimuthal 
Eq. [15] also corresponds to a parabolic potential well, p{6q) is also Gaussian with variance 
cr| = TqDq/2 which gives \im.dt^oo{[dOQ{dt)Y') = DqTq. Again, since Dq is measured from the 
slope of {[d9Q{dt)Y') ., one can obtain the time scale Tq = 6.9 s from the ratio of measured 
values. The transition between scalings occurs at dt = Tg, so this is the time scale over which 
^0 remains correlated. 

Allowing for variations in 5, for instance during cessations, the time scale corresponding 
to the damping term for Eq. [15] is Tg6o/6, which diverges when 6 = 0. However, such a large 
time scale does not exist for the LSC dynamics since cessations have a typical duration on 
the order of ts (see Sect. IVI CI) . Since 2ts is the correlation time for 6, the damping term 
of the azimuthal equation may retain some autocorrelation over this time scale. This is the 
longest time scale expected for the dynamics of 6*0, and in fact Fig. H] has a slightly lower 
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plateau for dt ^ ts. 

The data in Fig. Hlwere intended to test whether the fluctuations in 60 are diffusive, but 
the range of the sloped region is too short to do this with confidence. At best the diffusivity 
can be estimated based on the first two data points by assuming diffusive behavior. The 
plateau region of Fig. H] corresponds to a range where the dynamics can be considered 
strongly damped. For these time intervals, the acceleration term is negligible (^0 ~ 0) and 
the azimuthal Eq. [15] becomes 



^0 



(22) 



for the diffusion of ^o- If implies that fluctuations of 6q follow a diffusive scaling {d6o{dt)^) 
Dodt where the diffusivity is 



Dn 



6 



(23) 



Diffusive behavior for On was observed in Refs. 



21 



22 



23l | and was used in Ref. [2J] to study 



the long- 
in Refs. 



;erm dynamics of 9q due to Earth's Coriolis force. The variances {{dOY) reported 
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24f | do not saturate at a maximum value because 6*0 represented by Eq. [22] 



there are weak 



has no potential terms and thus it is unbounded. In the model of Ref. 
potential terms that would go on the right-hand-side of Eq. [22] due to Earth's Coriolis force 
which tend to align the flow in a preferred orientation and cause a net azimuthal rotation 
of the LSC. However, these terms are small compared to the fluctuation size, so they do 
not effectively bound 9q to a flnite range {2^. Since a weak Coriolis force is only relevant 
to long-term dynamics, a long-term approximation can be made by assuming 5 = where 
the predicted value of the diffusivity of Oq is Dq = r^Dg = 1.4 x 10~^ rad^/s. This can be 
compared with the value Dq = 1.22 x 10^'^ rad^/s reported in Ref. 24|. The good agreement 
shows that the Coriolis force model of Ref. 2^ is consistent with the strong-damping limit 
of Eq. [15] near the stable flxed point. Because the Coriolis-force terms are small they do not 
effect any of the short-time-scale dynamics of 6 or 60 and can be ignored in the present work. 
Thus the results of Ref. 2J] are completely consistent with and represent a long-time-scale 



limit of the current model. 



14 



C. The probability distributions 



1. The distribution p{5) 

The probability distribution p{S) of the amphtude S can be calculated from the steady- 
state solution of the Fokker-Planck equation (see, for instance, Ref. 4^) which represents a 
balance between advection and diffusion of probability: 

«.)^^^ (24) 
where 5d is the deterministic part of Eq. [TOl The solution to this differential equation is 

p{5) oc exp (25) 

where the potential Vs is given by Eq. [191 In the linear approximation, valid near the stable 
fixed point, the potential from Eq. [20]is parabolic, and then p{6) is Gaussian with variance 

a| = TsDs . (26) 

Figure [5] shows the probability distribution p{5) derived from experimental data at i? = 
1.1 X 10^° as open circles. The predictions of p{6) for both potentials Vs and are plotted 
as well as dotted and dashed lines respectively. For the predictions the values of Ds and 
obtained from Fig. [3] and of 6q derived from the experimental time series for 6{t) were used. 
The Gaussian shape of the peak and its variance are correctly predicted by the model. The 
good match near the peaks supports the validity of the linearized model of diffusion in a 
potential well near the minimum. While the nonlinear model correctly predicts a skewed 
distribution favoring small 6, the predicted p{6) does not match the experimental data in 
the tails of p{6) very well. 

For 6 < 0.55o, the measurements shown in Fig. [5] suggest an exponential dependence of 
p{6) on 5. A fit of 

p{5)=p{0)exp^ (27) 

with the free parameters B and p{0) to data for R = 1.1 x 10^° is shown in the figure. Since 
cessations occur due to large fiucutations in 6, the statistics corresponding to cessations 
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should be determined by the tail of p{S) rather than by p{S) near its peak. Because the 
model is based on the assumption of the existence of the LSC, it is perhaps not surprising 
that it fails when the LSC is weak. Near the end of this paper, in Sect. I VI It we suggest 
a possible expansion of the model that could reproduce the small-5 behavior of p{6) more 
accurately. 



2. The distribution p{Oq) 

The peak of the probability distribution of p{9q) can be estimated from the linearized 
azimuthal equation, Eq. [181 This equation is mathematically similar to the amplitude equa- 
tion Eq. [T71 and the steady-state solution of the Fokker-Planck equation leads to the similar 
result 

p(^o)cxexpf-^') (28) 
for 6*0 near the stable fixed point. Thus the variance oipiOo) is 

a| = ^ . (29) 

Figure [6] shows the experimental results for p{9q) for R = 1.1 x 10^° along with the model 
predictions for constant damping, as well as the result of a numerical simulation of the 
full model Eqs. [10] and [151 The linearized model accurately describes ^(6^0) as a Gaussian 
near the peak of the distribution, but with larger tails due to the variable damping term 
which allows large fluctuations in 6q when 5 is small. The simulation distribution does not 
match the tail of the experimental distribution very well, but does have an approximately 
exponential decay for large 9q like the experimental data. 



D. Power spectra 

To a limited extent one can test the assumption of white noise for the stochastic terms 
/^(t) and /g(t) by examining power spectra of 5 and 9q. These are shown in Fig. [3 for 
R = 1.1 X 10^°. Both power spectra are mostly flat for small frequencies, as expected for white 
noise, but they show a roUoff for large frequencies. The measured power spectra are a result 
not only of the stochastic terms but also of the response to them contained in the dynamical 
equations Eqs. [TOl and [151 Since the system spends most of its time near the stable fixed 
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point, the amplitude response can be calculated to a good approximation from the linearized 
equation Eq. [T71 For a single driving frequency /, the stochastic term can be represented 



by fsit) = y Ds/dtexp{i2TTft). Assuming responses of the form 6 = C exp{i27rft + z$), 
substituting these into Eq. \17\ and taking the magnitude of the complex solutions leads to 
a power Ps given by the Lorentzian function 

A fit of this function to the data very near the scaling crossover gives = 34 s and Dg = 
1.8 X 10~^ K^/s. Considering the approximations made in deriving the model, these values 
are in satisfactory agreement with those obtained from fits of {[d6{dt)]'^) (see Table H]). The 
power spectrum Ps has unexplained small features (which also appear in other data sets) that 
deviate from the expected function. In addition, the roUoff appears to have an exponent 
somewhat more negative than —2, indicating that there are more low-frequency and less 
high-freqeuncy fluctuations than expected. Nonetheless, the shape of Ps is at least roughly 
similar to the model prediction. 

Similarly, the response to the linearized azimuthal equation Eq. [15] is 

(31) 



A fit of Pg to the power spectrum of gives Tg = 5.4 s and Dg = 1.4 x 10~^ rad^/s^. This fit 
is excellent, and as was seen for Ps and {[ddldt)]"^) , the fit parameters are also in satisfactory 
agreement with those found from fits of {[d6o{dt)]'^) . This power spectrum is consistent with 
the model assumption of white noise which is filtered at high frequencies due to inertial 
damping. 

E. Correlation between 5 and |^o| 

So far all of the tests of the model have assumed the damping term of the azimuthal 
equation Eq. [15] to be constant in order to make the azimuthal equation linear. A need for 
a variable 6 to appear in the azimuthal equation was indicated already by measurements 

25^. 



that showed that 16*01 and S are negatively correlated with 5 slightly leading 16*01 [23, 



This is now explained qualitatively by the model; as the damping term increases with 5, 
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the magnitude of the azimuthal rotation rate decreases. Since ^ r^, changes faster 
than S. Thus the distribution of comes close to a temporary stationary state for a given 



6 after a time of the order of ta. This is the lead time measured in Ref. 



251 ] ■ It was found 



experimentally in Ref. 25|] that the delay time was about 6% of the turnover time T. Using 
Eq. [T6l and Re = 2L^/(Tz/) one can estimate roughly that Tq = T/4. The measured Tg is 
somewhat smaller, and it has a somewhat different dependence on R (see Fig. [8] below). In 
either case there is order-of-magnitude agreement with the measured lead time. Numerical 
simulations indicate that the peak of the correlation occurs at a time of the order of r^, but 
that it also depends on other model parameters. 

V. THE RAYLEIGH-NUMBER DEPENDENCE OF THE PARAMETERS 
A. The time scales ts and 

The time scales ts and Tq predicted by the model were inferred from the mean-square 
variable change over time, and from power spectra, as discussed above. The azimuthal 
time scale can also be estimated in the strong-damping limit from the ratio of measured 



diffusivities = ^JDq/Dq. Each of these is non-dimensionalized in the same way as the 
turnover time to obtain a Reynolds number. This gives -Rf = L"^ /{tsu) which is plotted in 
Fig. [Hk and -Rf = L"^ /{tqv) which is plotted in Fig [8b. The different methods of measuring 
the time scales all agree within about a factor of two. The model predictions of Eq. [12] 
and [in] are also shown in the figure. Fits of power laws to the Reynolds numbers obtained 
from variance measurements in the large sample gave i?f oc R^'^^ and -Rf oc R^'"^^. These 
exponents do not agree with the predicted exponents of about 1/4 from Eq. [12] and about 
1/2 from Eq. [TG] respectively. The Reynolds numbers for different samples with the same 
cr and F do not agree at the same i?, indicating that the L-scaling of the model prediction 
is incorrect. The Reynolds number R^ = 2LF'/{Tv) corresponding to the plume circulation 



period T was reported in Ref. 43| to be Re = 0.0345/?^^^ and is shown in Fig. [8^. It 
is seen to be very close to i?f for the large-sample auto-correlation measurements, which 
implies ts ~ T/2. The two predicted time scales were observed by several methods, and the 
predicted values are within an order-of-magnitude of the data, but both the L-scaling and 
the i?-scaling of the model disagree with the measurements. 
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In principle, the time scales could be affected by the stochastic terms from the Langevin 
equations if these terms are auto-correlated over some time interval. The time scales obained 
from the variances represent long-term dynamics. Thus they should be unaffected by the 
correlation times of any of the model terms. A power spectrum could be modified if the 
corresponding stochastic term has a non-zero correlation time, which corresponds to col- 
ored noise. Since the time scales measured from the correlation functions are close to those 
obtained by the mean-square variable change, the correlation times of fs{t) and fg{t) must 
be small compared to ts and respectively. The correct prediction of the peaks of prob- 
ability distributions based on the mean-square variable change provides support for using 
the timescales obtained by this method, so these measured timescales will be used as the 
experimental input for other model predictions. 



B. The amplitude 60 



Previous work revealed correlations between t 
near the side wall in the horizontal mid-plane [l7 



le velocity of the LSC and the temperature 
Equation [11] predicts the relationship 



between the mean temperature amplitude 60 and the Reynolds number. It can be rearranged 
to get 

A X ^ ^ 187r(i?,)i . (32) 

Figure [9] shows measurements of 60/ AT x R/a vs. Re over 2.5 decades of R and for 3.3 < 
a < 5.5, with data from both the medium and the large sample. Values of Re are based 



on the plume circulation period determined in Ref. [43j from measurements with the same 
apparatus as the current experiments. The Prandtl number a was varied by changing the 
mean temperature of the fluid. The solid line shows a power-law fit with the exponent 3/2 to 
the data with a free coefficient c defined by 5q/ NT x R/a = c{Re)^^'^. The fit yields c = 159, 
a factor of 2.8 larger than the prediced ISvr. This power law fits the data within better than 
20% over 2.5 decades of R, and thus strongly supports the predicted relationship between 
60 and Re- 
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C. The difFusivities 

The values of the non-dimensionahzed diffusivity Ds x (L^/AT^z/) are shown as a function 
of R in Fig. (TDl The non-dimensionahzation used there leads to a disagreement between 
data from the two samples for the same control parameters R and a, which is undesirable. 
A power law was fit to the large-sample data and yielded an exponent of -0.04. Since 
(AT)2 oc R\ this gives Dg oc R^-^^. 

The value of the non-dimensionalized diffusivity Dg x [L? /vY is shown as a function of 
R in Fig. [TTJ This non-dimensionalization is also seen to lead to disagreement between data 
from the two samples at the same R. A power law was fit to the large-sample data to obtain 



D. Non-dimensional parameters 

An alternate non-dimensionalization can be made by combining the three parameters 
from the equation for the temperature amplitude 5, Eq. [TOl into 7 = {Dsts)/6q, and the two 
parameters from the azimuthal equation, Eq. [151 into DqT^. These are shown in Figs. [121 
and [131 respectively, for various R in both samples. While there is still some disagreement 
for 7 between the two samples, it is smaller than in Fig. [10] and it is possible that this may 



be due to a non-Boussinesq effect 38|. These two dimensionless parameters have only a 
weak i?- dependence. A fit of a power law to D^r^ gives an exponent of 0.14 ± 0.01. In the 
large sample 7 is essentially constant, while in the medium sample and for i? ^ 4 x 10^ it 
varies as R^'^"^ as shown by the solid line in Fig. [121 These two non-dimensional parameters 
completely determine the parameter space of the linearized model. Taking into account 
the variable damping term requires the additional non-dimensional parameter t^/t^ which 
is proportional to R^-^^ in the large sample. Since all of these non-dimensional parameters 
vary only weakly with R, the large experimental range of R covers only a small region of 
the parameter space . It would be very interesting to learn how the parameters vary with 
the aspect ratio F of the sample. 
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VI. CESSATION RESULTS 



A. Empirical potential 

Cessations occur when 5 becomes small, so the approximations near the stable fixed 
point may not be adequate for describing them. It was seen in Fig. [5] that the model 
does not accurately predict the tails of p{S), so Eq. [10] will not accurately describe the 
statistics of cessations, even though cessations do occur for that dynamical equation at a 
rate which exceeds the experimental observations by only a factor of two or three. The 
problem can be seen directly in a comparison of the structure of the Langevin equation for 
6 with experimental measurements. It was observed that the average of 6 is independent of 
6 both before and after a cessation when 6 is sufficiently small, roughly when 6 < 6q/2. [23] 
For small 6, the model Eq. [10] predicts 5 oc 5, in disagreement with the data. In this section 
it will be shown that this inconsistency can be repaired by using an empirical potential 
inferred from the experimentally measured small-5 tails of p{S), grafted onto the parabolic 
potential of the model near 6 = 60. 

On its own, p{6) determines the product of the frequency and of the duration of cessations. 
These two values can be predicted separately for diffusion in this empirical potential well Vs^e 
by also using the model parameters that were obtained in Sect. IIVI The empirical potential 
is obtained by fitting Eq. [27] to the experimental results for p{6) for small 6. For small 6 the 
form of Eq. [27] and the Fokker-Planck result Eq. [25] then imply 

Vse = ^ + constant . (33) 

We note that this potential does not yield a fixed point at 5 = because its derivative 
is finite. Nonetheless it yields a well defined potential barrier that can be used to predict 
properties of cessations. The time-averaged 6 is directly related to the potential by the 
equation 6 = —dVs/d6. For small 6, Eq. [33] implies 6 = BDs/{26o) which is independend of 
6 in contrast to the model Eq. [TU] The potential barrier AVs^e = ^,e(0) — Vs^ei^o) can be 
expressed in terms of the measured p{6) using Eq. [25] to get 

^^J^ = lnp(5o) - Inp(O) = A Inp . (34) 
The values of B'j and 'jAlnp obtained from fits of p{6) = p{0) exp{B5/5Q) at various R 
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are shown in Fig. [TH The peak value p{So) = {2i:DsTs)~^I'^ is taken from the Gaussian 
approximation of p{5). Both values are seen to be nearly independent of R when scaled 
by 7. The inferred potential is shown in Fig. [2] for 5 < 0.5(5o using the typical measured 
values B = 8.6 and 7 = 0.069, and A Inp = 6.2 that are a good approximation for all of the 
large-sample data. The potential barrier is larger for the empirical potential then it is for 
the original model potential, and thus it yields a smaller cessation frequency. 

For the potential to be smooth and thus for Sd to be well defined in the corresponding 
Langevin equation, the small-(5 limit given by Vs^e and the potential Vs near 6 = 6q must 
match up at some intermediate value 6 = C6q. Based on the measured p{S), C ~ 0.5. These 
requirements fix the parameters B and A Inp to be i? = (1— C)/7 and A Inp = (1— C^) / (27). 
Thus the observation that i?7 and 7 A Inp are roughly constant is equivalent to C being 
roughly constant, and that p{5) has the same shape for different control parameters. Since 
we do not have data covering a large range of 7, it cannot be confirmed that both B and 
A Inp scale as I/7, but the fact that the significant decrease in 7 for small R shown in 
Fig. [12] does not appear in Fig. [14] is consistent with that conclusion over a small range of 7. 



B. The rate of cessations 



The equation for 5 does not contain an inertia term, and corresponds to pure diffusion in 
a potential well. Thus, if the root-mean-square amplitude of this diffusive motion is small 
compared to its mean value ^o, then successive cessations will follow Poissonian statistics. 
This will be the case when the diffusivity Ds is small compared to the depth of the well, 
and this condition is satisfied for the physically relevant parameter values. 

The experimental results for the time-averaged frequencies of cessations tUc, measured in 
events per unit time, are not very accurate because cessations occurred only about once 
or twice per day. Nonetheless a comparison between experiment and the model is useful. 
The rate of cessations can be calculated using the model of diffusion across a potential 
barrier, which is analagous to the Arrhenius-Kramers problem 42]]. Stochastic fluctuations 
with strength Ds drive the amplitude 5 in a potential well. A cessation occurs when 5 
fluctuates from the bottom of the potential well aX 5 = 5q to the top at 5 = 0. Thus, for 
a cessation to occur, fluctuations must overcome a barrier = V{0) — V{6q). First we 
shall use the empirical potential barrier AVs^e- It can be obtained from Eq. [M] and the 
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data in Figs. [12] and [TH One sees that the ratio of the potential barrier to the diffusivity 
given by Ds is around 6 or larger. This means that the Arrhenius-Kramers equation 

for the large-barrier approximation 4^ can be applied here. In that approximation the 
rate of escape is proportional to exp(— 2AVi5^e/D5). The prefactor depends on expansions of 
integrals around the peak and minimum of p{5) with the result depending on the curvature 



of the potential at these points. The full solution is Uc = {B / ts) y / {2^) exp{—2AVs,e/Ds)- 
While the exponential dependence always has the same form in the large-barrier limit, the 
proportionality to B^/^ is a result of the shape of the potential near 6 = 0. This prediction 
for ujc is shown in Fig. [15] as triangles, along with experimental data from Ref . 23|] which are 
given as circles. The prediction is seen to be in agreement with the experimental data within 
about a factor of 2. The cessation rate is nearly independent of R for the large sample, and 
decreases for the medium sample with decreasing R. This plot roughly follows the same 
trend as 7, see Fig. [121 confirming that 7 is the most relevant parameter for determining 
the rate of cessations. This it must be since it is the only dimensionless parameter in the 
relevant Langevin equation. 

For comparison, the prediction in the large-barrier limit for the model potential Eq. [T9] is 
uJc = exp{—2AVs/Ds) / {27rTs), but the large-barrier limit is not as good an approximation for 
this model in the experimental parameter range because AV^ is too small (see the potential 
barriers in Fig. [2|). For instance, for R = 1.1 x 10^^ it yields about 14 cessations per day, 
which is an order of magnitude larger than the experiment. However, when the complete 
model, Eqs. [TO] and [H] is integrated numerically for the same R instead of using the large- 
barrier analytic expression, one finds about 3.8 cessations per day [33], which is only a factor 
of two larger than the experiment. It turns out that the prediction based on the empirical 
p{6) and the prediction based on the parabolic potential of the linearized model give about 
the same result for Uc, but the functional dependence on 7 differs between the two due to 
the different shape of the potentials around 6 = 0. 



C. The duration of cessations 



For 6 ^ 0.55o it was observed experimentally 23]] that during cessations the average over 
all events of the magnitude {\6\) of the rate of change of 6 was independent of 6. Thus the 
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amplitude drop and recovery are symmetric and follow the equation 

S{t)=6m + m\t\ (35) 

where t is the time elapsed since the cessation. We calculate the average duration of a 
cessation to be (At) = {60 — 26m)/ (|^|)- Here 6m = 0.0955o is the average minimum measured 
amplitude during cessations, which is just slightly larger than zero in any experimental 
measurement. The value of {6 (6)) for the rise can be calculated from the experimentally 
obtained small-5 tail of p{6) to be (6) = —dVs^e/d6 = BDs/{26q) using the empirical Vs^e from 
Eq. [33l This (6) is independent of 6 in agreement with experiment [23] . The average duration 
of cessations is thus predicted to be (At) ^ {60 — 26m) / [B Ds / {26o)]. The average measured 
cessation duration, as well as the predicted value based on the empirical potential, are shown 
in Fig. [in] for various R. Both data and the prediction of At are roughly proportional to r^, 
but the prediction is an overestimate by about a factor of 2. While the amplitude behavior 
during cessations is inconsistent with the proposed model Eq. in the sense that it is 
responsible for the tails of p{6), using the tails of the measured p{6) as experimental input 
into the diffusion model allowed the prediction that (6) during cessations is constant and 
approximately proportional to 60 /rs. 



D. Angular change during cessations 

Within experimental resolution the net angular change A6 during cessations was found 
to have a uniform probability distribution p{A6) 25|]. Numerical simulations based on 
the model equations Eqs. [TU] and dSl confirmed this within their resolution 33|| . A viable 
but ad hoc explanation of this result was that, once the LSC ceases, there is no memory 
of its original orientation and the re-organization of the new circulation will occur in an 
arbitrary new orientation. Here we offer an alternative view of this phenomenon. We 
present arguments showing that a near-uniform distribution can be due to large azimuthal 
fluctuations that occur when the inertial damping term in Eq. [15] becomes small. This in 
occurs only when 6 becomes small during cessations. 

Equation [28] predicts that, over a given time period 6t, the orientation diffuses, yielding a 
Gaussian distributed p{A6). The typical orientation change due to diffusion near the fixed- 



point amplitude 60 is d6rms = VDeAt, which with Eq. [15] yields dOrms ~ DnAt. For 
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the average cessation period {At) ~ 1.34T [23| one then has dOrms ~ 0.3 rad. Because the 
inertial damping decreases as 6 decreases during a cessation (see Eg. [T5l) . the angular change 
is typically larger during cessations then it is when 6 is near 6^ or larger. In the strong- 
damping approximation, the diffusivity of 6q is given by Eq. [23] and is larger at smaller 6, 
and the time dependence of 6(t) is given by Eq. [33 The mean-square change in amplitude, 
when integrated over the duration of a cessation, gives 



/•At/2 

{[deo{At)]^) = J_^^ ^ D0[6{t)]dt = 2D^rl6oAt/6^ . (36) 



rAt/2 

'-At/2 

Thus the typical angular change during cessations a^e = \/ {[dOoiAt)]"^) ^ lA rad at i? = 
1.1 X lO^''. This is much larger than the angular change for constant damping with 6 near 
6o that was illustrated above. 

There is an additional contribution during cessations when the flow reverses - that is 
when 6 becomes less than zero. This can occur due to a continuous meandering of 6 below 
zero, but in our analysis the absolute value \S\ is taken and ivr added alternatively to 6q 
whenever a reversal occurs. For an odd number of reversals, there is a net contribution of 
vr to A6. The percentage of cessations with reversals depends on the potential barrier to 
cessations and how the threshold for cessations is defined. This would be 50% if a cessation 
was counted only if 6 crossed zero, but because the LSC cannot be resolved experimentally 
when 5 ^ 0, a cessation is counted when 6 drops below 0.155o- The chance of getting 
an odd number of reversals in a cessation can be estimated in the large-barrier limit to 
he A = 0.5exp[2[\/5,e(0) - Vs,e{0.l55o)]/ Dg] ^ 0.5 exp(-0.155) ^ 0.14 at R = 1.1 x 10^°. 
While this is at best an order-of-magnitude estimate, numerical simulations indicate that 
there are an odd number of reversals for a fraction A = 0.32 of cessations at this value of 
R. The combination of diffusion and reversals results in a double-peaked distribution at 
and vr for small azimuthal fluctuations, but with larger fluctuations, both peaks spread 
out. This predicted distribution is reduced to the range to vr by the transformation 
AOred = TT — Itt — \A6 mod 27r|| so that A6red is the smaller of the choices of angular 
change in either direction during the cessation. This transformation is made because of 
the non-uniqueness of 6o in which a change of and A^^ ± 27r are indistinguishable unless 
the orientation can be smoothly traced in time, which cannot practically be done during 
cessations. The distribution can be expressed as 
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p{Mred) OC ^ (1 - A) exp 



n=—oo 
oo 



+ ^ Aexp 



(37) 



If both peaks are assumed to have the same height, the net distribution is within 10% of 
uniform if a as > 1-35 rad. In the limiting case where A = (no reversals), for p{A6) to be 
within 10% of uniform would require a^e > 2.70 rad. For stronger diffusion, the absolute 
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Ref. 



ar change can increase, but p{A6red) remains nearly uniform. Experimental data from 
23| for 10^ < i? < 1.1 X 10^'^ are plotted in Fig. [TTl along with the predictions of Eq. [33 
for i? = 1.1 X 10^° and R = 10^ using A = 0.32. The model result is consistent with the 
data without any fit parameters. The predictions for both ends of the R range, as well as 
the uniform distribution, are consistent with the data because the experimental probability 
distribution has fairly large error bars since cessations are so rare that only a few hundred 
have been measured. 



Experiments reported in Ref. [45| with 1 = 1/2 for 1.5 x 10^° < < 7.2 x 10^° found 
p{A6) peaked at and vr with p{'k/2) ^ 0.3p(0). Although the effect of changing F on the 
model parameters is unknown, a single LSC roll was observed, so the model should apply 
to those experiments with different parameters. The measured p{A6) is consistent with the 
predicted functional form of the sum of Gaussians, one centered at each integer multiple of 



TT. 



VII. AN EXPANDED MODEL FOR SMALL 6 

Measurements of p{6) and 6 for small 6 indicate that the Langevin equation for 6 should 
be dominated by a constant term for 6 < 0.56o, in contrast to the model equation Eq. [10] 
which has 5 ~ 5 when 6 is small. In this section we outline how the model might be expanded 
so as to account for the small-5 behavior, but leave the details of this expansion to future 
work. 

The expansion consists of the inclusion of the heat-transport equation 

T + M ■ VT = kV^T (38) 
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to describe the top and bottom boundary-layer heat-conduction. When the LSC is weak, 
6 and U are small and the advective term of Eq. [38] can be neglected. The diffusion term 
contributes only in the thermal boundary layers and can be estimated as kV^T ^ AT/{2Q, 
where / is the width of one boundary layer. It is given approximately by / = L/{2M) [37,] . 
This term contributes in the fractional volume 21 /L. Using Eq. [31 the volume averaged 
temperature change is T = 4(5/(37r). Combining these terms gives 5 = 37tk,AT^//{2L'^) ^ 
0.03 K/s for R = 1.1 x 10^°. This is a constant, i.e. independent of 6, as required by the data 
for p{S). It is within an order-of-magnitude of the experimental value S ~ I.SASq/T ^ 0.007 
K/s reported in Ref. 23|. Using the result Af oc R^^^ for the large sample 46|, the prediction 
gives 5 oc R'^/^. This is close to the measured 5 (x 5o/T (x R^^^ reported in Ref. |23|. These 
values are consistent with the conclusion that the dominating driving term for small 6 is 
due to thermal diffusion across the boundary layers. At least for small 6, this replaces the 
need to use the momentum equation and the assumption 6 oc U. 

A more complete model that applies for all ranges of S might consider the LSC velocity U 
separately from 6. This separation would be more realistic since 6 and U fluctuate separately, 
even though the model produced cessations and rotations without taking this into account. 
Such a model would likely include the momentum Eqs. [7|and[T4|as well as the heat transport 
equation, for a total of three differential equations for the three parameters U, 6, and 6. The 
advective term of the heat-transport equation couples it to U, while the buoyancy term of 
Eq. [H for the LSC velocity couples it to 6. The fact that only one equation was necessary 
to describe the LSC strength near the stable fixed point suggests that the time scale for the 
coupling of 6 and U is short compared to rs, so that the dynamics of 6 and U can to lowest 
order be described by a single equation. 



VIII. SUMMARY AND CONCLUSIONS 



We demonstrated that many aspects of the dynamics of the LSC can be described by 
two coupled stochastic ODEs that are motivated by the Navier-Stokes equations. One of 
the equations is for the amplitude of the azimuthal temperature variation 6 that is induced 
by the circulation, and the other is for the azimuthal orientation Oq of the near-vertical 
LSC circulation plane. The 5-equation represents the balance between the driving due to 
buoyancy and the dissipation due to the drag across viscous boundary layers near the walls. 
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The 9o equation is coupled to the 5-equation by an advective term that represents the 
angular momentum of the circulation. Both equations are driven by a Gaussian white noise 
term that represents the action of the small-scale turbulent fluctuations on the large-scale 
circulation. 

The original definition of reorientations as a relatively large and fast angular change was 
cumbersome because there was no clear distinction between small fluctuations and large 



reorientations 



251]. Now all of the azimuthal dynamics - including meandering, larger and 



faster rotations, and the large azimuthal changes associated with cessations - are simply 
described as diffusive fluctuations in a potential well Vg given by Eq. [2T] and with a curvature 
that depends on 6. Larger rotations tend to occur when 6, and thus the well curvature, is 
smaller, and the large angular change during cessations (with 6 near zero) occurs because 
the potential well becomes nearly flat so that the azimuthal motion is almost unconfined. 
Reversals of the LSC with ~ tt are not distinct events in this description. Crossings of 
the potential maximum at 5 = tend not to be reversals because this is when the azimuthal 
fluctuations are at their largest. 

Previously, we had suggested that the uniform distrubution of p{A9) for cessations occurs 



because the LSC loses any memory of its orientation during cessations [23|. This is the 
expected result if the LSC completely breaks up. In Sect. IVIDI we showed that a good 
approximation to this distribution can occur even if 6 remains finite during cessations, 
provided 6 becomes small enough so that the azimuthal diffusion is large over the durations 
of cessations. 

By studying power spectra and probability distributions, it was found that the stochastic 
driving terms can to a good approximation be described by Gaussian white noise. The model 
contains two time scales: Tq and ts. These were measured indirectly by several methods to 
be within a factor of two or so of the prediction; however, the predicted dependence on 
the Rayleigh number R and the sample height L differed somewhat from the data. The 
prediction for the dependence of the mean temperature amplitude 6o of the LSC on R 
agreed well with measurements over the large range of R and the small range of the Prandtl 
number a that were explored. 

The model accurately describes the dynamics and much of the /^-dependence using the 
potential near the minimum where 5 ~ 5o and predicts the existence of cessations near 6 = 0. 
However, it fails to predict quantitatively the behavior far from the potential minimum. 
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including the tails of p{S) and details about cessations. Using the measured p{S) as empirical 
experimental input to the model, the scaling behavior of the frequency, angular change, and 
amplitude behavior during cessations could be reproduced quite well. This supports using 
a model of diffusion in a potential well for the LSC even though the originally proposed 
model equations must be modified for small 6. In the preceding section we suggest that the 
modification may consist of the inclusion of the driving by the diffusive heat transport across 
the top and bottom thermal boundary layers which becomes important when the driving 
due to 6 becomes small. 

In agreement with experiment, the model predicts a rate of cessations that is roughly 
uniform at about one to two per day over a wide range of Rayleigh numbers, from 3 x 10^ < 



R<10 



11 



231]. On the other hand, the model of diffusion over a potential barrier suggests a 



very sensitive dependence of the rate of cessations on the model parameters, at least when 
the potential barrier is not small compared to the diffusivity. The lack of a dependence 
on R can be understood because the barrier ratio 2V/Ds oc is nearly independent of 
R. It is conspicious that the two dimensionless parameters 7 and DgT^ that determine the 
dynamics of 6 and 6q respectively are both nearly constant over the 2.5 decades of R that 
were studied. 

The only known dynamical behavior of the LSC that is conspicuously missing from the 
model is a description of the twisting oscillation of the LSC [o], [2^. There is no signal of 
the twisting oscillation in measurements at the mid-height of the sample, so it seems that 
the oscillations are independent of the azimuthal diffusion and re-orientations described by 
the model. The Nusselt number, or non-dimensional heat flux, is another important aspect 
of RBC not represented in the model. The Nusselt number is understood to be controlled 
by the thermal boundary layers and not the LSC, so its absence from a model of the LSC 
which does not include these boundary layers is not surprising. 

Here the LSC was studied in a container with cylindrical symmetry. A previous paper 
described measurernents and a model of the symmetry-breaking effect of Earth's Coriolis 
force on the LSC [24]. That model is consistent with the strong-damping limit of the 
current model. A later project will consider the effects of various asymmetries of the system 
on the LSC in which perturbative terms can be added to the model equations Eqs. [10] and 

m 
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{[dx(dt)]'^) Px{^) prediction 
Ts (s) 47 34 85 

Tg (s) 6.9 5.4 13 

Ds (KVs) 6.4 X 10-5 1.8 X 10"^ 
(radVs^) 2.9 x 10"^ 1.4 x lO""^ 

TABLE I: The time scales Tg and r^, together with the model predictions Eqs.[T2]and[T6l and the 
diffusivities Ds and Dg, obtained from the mean-square change {[dx{dt)]'^) of the variables over 
time and from power spectra Px{uj) for R = 1.1 x lO^''. 
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FIG. 1: A diagram of the LSC showing the coordinate 0, maximum of the velocity profile U, and 
viscous boundary-layer width A (not to scale). 

FIG. 2: Solid line: The potential well V = [-d'^/2 + 26^/'^ /{5^/5^)]/Ts in which 6 meanders dif- 
fusively. Dashed line: a parabolic potential well = V{So) + {S — 5q)'^ /{Ats). Dotted line: the 
potential for 5 < O.5(5o inferred from the measured p{5). is the predicted potential barrier for 
cessations. 

FIG. 3: The mean-square change in amplitude {{dS)"^) as a function of the time interval dt for 
i? = 1.1 X 10"'^''. Solid line: a fit of {{d5)'^) = D^dt to the data for intermediate dt gives the 
diffusivity Ds- Dotted line: a constant 2Dsts fit to data with large dt, together with Ds, yields ts- 

FIG. 4: Dots: the mean-square change in azimuthal rotation rate {{dO)^) as a function of the time 
interval dt for i? = 1.1 x 10^*^. Open circles: modified time series using only azimuthal steps dO^ 
when 0.9(5o < 5 < \.\5q. Solid line: a fit of {dO^) = Dgdt to the open circles for dt < 6 s. Dotted 
line: a constant D^Tq fit to data with large dt. Together with Dg it yields . 

FIG. 5: The probability distribution of the amplitude p{5) for i? = 1.1 x lO^'^. Dashed line: 
Gaussian fit to data. Dotted line: model prediction from Fokker-Planck equation. Solid line: an 
exponential fit of p{d) = p{0) exp{BS/So) to data for 5 < O-GSq. 

FIG. 6: The probability distribution of the azimuthal rotation rate p(|^o|) over a single time step 
dt 2.5 s for i? = 1.1 X 10-'^'^. Circles: experimental data. Triangles: simulation data. Dotted 
line: the Gaussian distribution predicted for a constant damping term with 5 = Sq. Solid line: 
exponential fits to the tails of the distributions. 

FIG. 7: The power spectra Pg (dotted line) and (thin solid line) derived from the experiment. 
Dashed line: a fit of a Lorentzian near the crossover of Pg. Thick solid line: a fit of a Lorentzian 
to Pa. 
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FIG. 8: (a) The Reynolds number = L'^/{uts). (b) The Reynolds number ijf = L'^/{uTg). Each 
are measured by several methods for a wide range of R. Circles: from the mean-square variable 
change over time. Triangles: from power spectra. Squares: from the strong-damping approximation 



of Tq = ^Dq/Dq. Solid symbols: medium sample. Open symbols: large sample. Solid lines: 
prediction from model Eqns. [I2] and [T6l Dashed line: the Reynolds number corresponding to the 



turnover Re = 0.0345i?^/^ from Ref. 



43]. Dotted lines: power law fits to the open circles. 



FIG. 9: The measured value of (5) /AT x R/a as a function of Re- Solid circles: medium sample, 
a = 4.4. Open circles: large sample; a = 4.4. Up-pointing triangles: medium sample, a = 5.5. 
Down-pointing triangles: medium sample, a = 3.3. Solid line: a fit of the predicted power law 
(5) /AT X R/a = cRI^^ to ah of the data. 

FIG. 10: The non-dimensionalized diffusivity Ds x (L^/AT^z^) for various R. Solid circles: medium 
sample. Open circles: large sample. Solid line: a power-law fit to the large sample data. 

FIG. 11: The non-dimensionalized diffusivity X (L'^/i/)^ as a function of i?. Solid circles: medium 
sample. Open circles: large sample. Solid line: a power-law fit of the large sample data. 

FIG. 12: The non-dimensionalized diffusivity 7 = Dsts/6q for various R. Solid circles: medium 
sample. Open circles: large sample. Solid line: a power law with an exponent of 0.32. 

FIG. 13: The non-dimensionalized diffusivity D function of R. Solid circles: medium 

sample. Open circles: large sample. Solid line: power law fit to the data. 

FIG. 14: The coefficients obtained from a fit p{6) = p{0) exp^BS/do) to the measured p{5) for 
various R. Circles: Bj. Triangles: 7Alnp = 7 x [lnp((5o) — Inp(O)] = 2jV/Ds. Solid symbols: 
medium sample. Open symbols: large sample. 



23( 1 for various R. 



FIG. 15: Black circles: The measured frequency of cessations from Ref. 
Triangles: prediction using the measured tail of p{6) in the large-barrier limit for diffusion in a 
potential well. Solid symbols: medium sample. Open symbols: large sample. 
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FIG. 16: The average duration of cessations {At)/Ts for various R. The duration is calculated as the 
time that 6 remains below the threshold 5o/2. Circles: experimental data. Triangles: predictions 
based on measurements of the tails oi p{6). Solid symbols: medium sample. Open symbols: large 
sample. 



FIG. 17: The probability distribution p(A^red) of the orientation change during cessations, reduced 
to the range 0...7r by the transformation AOred = tt — \tt — \ A9 mod 27r||. Solid circles: data from 



23 1 . Solid line: uniform distribution. Dotted line; prediction from Eq. [37] for R = 10^. Dashed 
line: prediction from Eq. [37]for ii = 1.1 x 10^'^. 
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Fig. 1, Eric Brown, Physics of Fluids. 
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Fig. 12, Eric Brown, Physics of Fluids. 
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